/*
 * lib.h
 *
 *  Created on: 2010.05.10
 *      Author: Tadas
 */

#ifndef LIB_H_
#define LIB_H_

#define EPS 1.0E-6
#define JMAX 30
#define K 5

using namespace std;

static double tcount = 0;

double trapecY(double(*func)(double, double), double(*fna)(double),
		double(*fnv)(double), double x, int n) {
	double y, tnm, sum, del, a = fna(x), b = fnv(x);
	int it, j;
	if (a == b) {
		return 0.0;
	}
	static double s;

	if (n == 1) {
		return (s = 0.5 * (b - a) * (func(x, a) + func(x, b)));
	} 
	else {
		for(it = 1, j = 1; j < n -1; j++){
			it <<= 1;
		}
		tnm = it;
		del = (b - a)/tnm;
		y = a + 0.5 * del;
		
		for(sum = 0.0, j = 1; j <= it; j++, y += del){
			sum += func(x, y);
		}
		
		s = 0.5 * (s+(b-a)*sum/tnm);
		
	}
	
	return s;
}

double sumTrapecY(double(*func)(double, double), double(*fna)(double),
		double(*fnv)(double), double x) {
	int j;
	    double SNew, SOld = -1.0E-31;

	    for( j = 1; j < JMAX; j++){
	    	SNew = trapecY(func, fna, fnv, x, j);

	        if(fabs(SOld - SNew) < EPS*fabs(SOld)){
	            return SNew;
	        }

	        if(SNew == 0.0 and SOld == 0.0 and j > 6){
	            return 0.0;
	            // return SNew;
	        }

	        SOld = SNew;
	    }

	    cout << "Per daug iteraciju; rezultato nera" << endl;
	    return 0.0;
}

double trapecX(double(*func)(double, double), double(*fna)(double),
		double(*fnv)(double), double a, double b, int n) {
	tcount++;
	double x, tnm, sum, del;

	static double s;
	int j, i, it;

	//int it;

	if (n == 1) {
		return (s = 0.5 * (b - a) * (sumTrapecY(func, fna, fnv, a))
				+ sumTrapecY(func, fna, fnv, b));
	}

	for (it = 1, i = 1; i < n - 1 and i < JMAX; i++)
		it <<= 1;
	tnm = it;
	del = (b - a) / tnm;

	x = (a + 0.5 * del);

	for (sum = 0.0, j = 1; j++ <= it; x += del) {
		sum += sumTrapecY(func, fna, fnv, x);
	}

	return (s = 0.5 * (s + sum * del));

}
/// Kazkur cia da**ja kodo truksta
double sumTrapecX(double(*func)(double, double), double(*fna)(double),
		double(*fnv)(double), double a, double b) {
	int j;
	double SNew, SOld = -1.0E-31;

	for (j = 1; j < JMAX; j++) {
		SNew = trapecX(func, fna, fnv, a, b, j);

		if (fabs(SOld - SNew) < EPS * fabs(SOld)) {
			return SNew;
		}

		if (SNew == 0.0 and SOld == 0.0 and j > 6) {
			return 0.0;
			// return SNew;
		}

		SOld = SNew;
	}

	cout << "Per daug iteraciju; rezultato nera" << endl;
	return 0.0;

}

#endif /* LIB_H_ */
